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ABSTRACT 


The  finite  difference  (FD)  method  yields  the  solution  to  a  discretized  version  of  the  full  acoustic  wave  equation  for 
arbitrarily  complex  media.  It  is  a  full  spectrum  approach  and  is  thus  reliable  at  all  angles  of  propagation,  including 
backscatter.  This  offers  an  advantage  over  other  standard  propagation  methods  in  wide  use,  as  it  allows  for  accurate 
computation  of  acoustic  energy  levels  in  the  case  where  significant  scattering  can  occur  near  the  source,  such  as  may 
happen  for  an  explosion  near  the  surface,  or  underground.  This  fits  in  with  nuclear  monitoring  goals,  in  that  it  allows 
for  an  improved  understanding  of  the  generation  and  propagation  of  infrasound  energy  from  arbitrary  sources, 
including  underground  and  near-surface  explosions. 

Two  types  of  FD  methods  of  solving  the  acoustic  wave  equation  are  presented  in  this  paper.  The  first  is  a  finite 
difference  frequency  domain  (FDFD)  method,  applied  in  cylindrical  coordinates  to  simulate  the  effects  of  a  point 
source  in  an  azimuthally  symmetric  medium.  The  second  is  a  finite  difference  time  domain  (FDTD)  approach 
including  the  effects  of  both  gravity  and  wind,  applied  in  two-dimensional  Cartesian  coordinates.  In  this  paper 
equations  are  developed  for  the  FDTD  approach  where  both  wind  and  gravity  effects  are  considered. 

It  is  shown  that  the  FD  approach  can  be  used  to  solve  for  sound  intensities  in  arbitrarily  complex  models  that  may 
include  high  material  contrasts  and  arbitrary  topography.  In  this  paper,  results  of  FDTD  and  FDFD  approaches  are 
compared  for  the  case  of  a  shallow  underground  source,  for  a  boundary  with  significant  topography.  The  effects  of 
wind  and  gravity  on  the  solution  are  examined. 
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OBJECTIVES 


An  FDTD  method  of  numerically  synthesizing  infrasound  energy  in  a  realistic  environment  is  sought.  The  FDTD 
method  yields  the  solution  to  a  discretized  version  of  the  full  acoustic  wave  equation  for  arbitrarily  complex  media. 
It  is  a  full  spectrum  approach  and  is  thus  reliable  at  all  angles  of  propagation,  including  backscatter.  This  offers  an 
advantage  over  other  standard  propagation  methods  in  wide  use,  as  it  allows  for  accurate  computation  in  cases 
where  significant  scattering  can  occur  near  the  source,  such  as  may  happen  for  an  explosion  near  the  surface  or 
underground.  The  effects  of  wind  or  gravity  on  infrasound  propagation  may  also  be  incorporated  into  the  infrasound 
propagation  problem  with  relative  ease  using  finite  difference  techniques.  This  fits  in  with  nuclear  monitoring  goals, 
in  that  successful  completion  of  this  project  will  allow  for  an  improved  understanding  of  the  generation  and 
propagation  of  infrasound  energy  in  arbitrarily  complex  environments. 

In  this  paper,  both  an  FDFD  approach  and  an  FDTD  approach  are  used  to  solve  for  the  infrasound  signals  generated 
by  an  underground  source.  A  method  of  incorporating  wind  and  gravity  are  outlined  for  the  FDFD  approach. 

RESEARCH  ACCOMPT  JSHED 


An  FDTD  method  of  numerically  acousto-gravity  waves  in  a  windy  environment  has  been  developed.  The  equations 
governing  infrasound  propagation  are  derived  below  for  the  time  domain;  their  incorporation  into  an  FDTD  method 
are  explained  in  the  next  subsection.  Examples  are  given  in  the  final  sub-section  and  compared  to  the  frequency 
domain  method. 

Low  frequency  wave  equations  for  a  fluid  in  motion 

In  the  absence  of  viscosity,  the  equations  governing  propagation  of  sound  in  the  atmosphere  are  the  conservation  of 
momentum. 


conservation  of  mass 


and  the  ecination  of  state 


DV 

,_  =  _VP  +  F, 


Ri 

Dt 


-h  pV  •  V  =  0, 


2£p 

Dt  Dt  ’ 


(1) 


(2) 


(>•5) 


(Gill,  1!)82:  Ostashev  et  ai,  2005).  These  equations  relate  the  velocity  V,  the  pressure  P,  and  the  density 
p.  The  external  forces  acting  upon  the  medium  are  denoted  by  F.  and  C  is  related  to  the  adiabatic  sound 
speed.  At  low  freqiiencies.  the  gravitational  force  F  =  — pg  must  be  included,  where 


g  =  [0.0, 9.8]  m/.fl 


(4) 


indicates  gravitational  acceleration:  the  negative  sign  on  the  force  indicates  that  a  downward  force  acts  ujjon 
a  irositive  density  fluctnation  caused  by  the  propagating  sound  wave. 

The  convective  derivative  (also  known  as  the  Lagrangian  derivative)  ^  is  defined  by 


_D 


Dt 


+  (V  •  V). 


(r.) 


The  derivative  on  the  left  represents  the  change  with  time  in  a  reference  frame  moving  with  the  fluid.  The  first  term 
on  the  right  side  of  Equation  (5)  represents  the  change  at  a  point  fixed  in  space.  The  second  term  represents  the 
change  as  the  observer  moves  with  the  fluid  at  the  velocity  V,  and  is  called  the  advective  term.  Generally,  quantities 
are  expressed  in  terms  of  a  fixed  point  in  space  in  order  to  compare  computational  results  with  observations  made  at 
stationary  sensors. 


The  propagation  of  sound  waves  in  the  atmosphere  introduce  fluctuations  in  the  pressure,  density,  and  velocity 
fields.  The  standard  procedure  in  solving  Equations  ~(l)-(3)  is  to  consider  a  solution  of  the  form 
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p  =  Po  +  Ps;  p=  p  o  +  p  s;  V  =  w  +  v;  =c^  +(c^)';  (6) 

where  Po,  p  o,  and  w  are  the  ambient  solutions  in  the  absence  of  fluctuations  Ps,  p  s,  and  v  due  to  the  propagating 
sound  wave.  The  adiabatic  sound  speed  is  given  by  c,  and  (c^)'  is  the  perturbation  in  the  squared  sound  speed  caused 
by  the  passage  of  the  sound  wave  (Ostashev  et.al.,  2005).  In  what  follows,  w  denotes  wind  velocity  profde,  and 
V  denotes  the  acoustic  particle  velocity  associated  with  the  sound  wave.  Waveforms  are  derived  by  computing  the 
pressure  perturbations,  Ps,  as  a  function  of  time. 


To  obtain  first  order  linear  equations  feasible  for  implementation  in  an  FDTD  method,  it  is  assumed  that  the  ambient 
field  is  in  equilibrium,  and  that  the  sound  wave  pressure  is  very  small  in  comparison  with  the  ambient  pressure.  The 
former  assumption  implies  that  Eqs.  (1-3)  are  satisfied  in  the  absence  of  a  propagating  sound  wave  so  that  zeroth 
order  terms  cancel  when  Eqs.  (6)  are  substituted  into  Equations  (l)-(3).  The  latter  assumption  means  that  second  and 
higher  order  terms  cancel.  Retaining  only  linear  terms.  Equation  (1)  becomes 


/  \  I  \ 

Po^^-^-bv  - Vw-bw  - Vvj  t  w  •  Vwj  =  -Vp„  -p.sg. 

where  the  external  gravitational  force  has  been  iiic'lnded. 


'()w 


(7) 


Under  the  assumption  that  winds  are  .steady,  i.c...  Ow/Ot  =  0,  and  there  is  no  variation  along  the 
direction  of  the  wind,  i.e.,  w  •  Vw  =  0.  the  ecpiation  governing  acoustic  particle  velocity  reduces  to 


/ilv 

\m 


)  v  •  Vw  I  w  ■ 


-Vp.,  -  p,,g. 


(8) 


Combining  Eqs.  (2)  and  (d),  and  inserting  Ecjs.  ((i)  yields  the  following  equation  for  the  pressure 
fluctuations  due  to  the  propagating  sound  wave: 

^  -t-  w  •  Vp,,  +  v  •  Vp,,  =  -pc^V  ■  v  -  (^Po{c'f  +  p„c^J  V  ■  w  (!)) 

The  hydrostatic  ecjiiation  Vpo  =  —pog  holds  under  the  assumption  of  steady  state  and  range-invariant  winds. 
Furthermore,  Ostashev  ef  al.(200r))  states  that  the  divergence  teriii  V  •  w  may  be  ignored  to  order 
under  this  a.ssnniption  Eq.  (!))  becomes 


-b  w  •  Vp.,  =  p„v  ■  g  -  poc  V  ■  V, 


where  the  gravity  vector  is  as  defined  in  Eq.  (4). 


(10) 


The  equation  for  flnctiiations  in  density  may  be  derived  from  Eq.  {,'{),  the  ecpiation  of  .state,  as 


Ops 

c)i 


t  V  ■  Vpo  +  w  •  Vp„  = 


1  fOp„ 


(II) 


under  the  assumptions  that  the  ambient  pressure  and  density  values  do  not  vary  with  time,  and  that  (c')^w- 
Vpo  =  0,  which  is  satisfied  by  horizontal  winds,  and  densities  that  vary  weakly  with  range.  This  ecpiation 
for  density  variations  is  similar  in  form  to  Ecj.  (i.l4.d  of  Gill  (1982),  with  the  exception  that  horizontal  winds 
w  are  included  in  Eep  11.  Making  use  of  the  hydrostatic  ecpiation  Vpo  =  — P.9,  the  above  ecpiation  may  be 
re-written  as 


dps 

(It 


v^poN'^lg 


dPs  s 
(It  > 


(12) 


where  (j  =  is  the  vertical  component  of  gravity,  v-  is  the  vertical  particle  velocity,  d/dt  represents 

the  time  derivative  in  a  reference  frame  moving  with  the  wind,  and  the  Brmit-Vaisala  frequency  N  -  also 
called  the  buoyancy  frecpiency  -  is  defined  as 


A' ^  =  -(j{Pi)  ^dpo/dz  +  cy/c^). 


(Id) 
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For  a  stable  medium,  N  is  a  positive  real  number  equal  to  the  frequency  of  a  parcel  in  vertical  motion.  Equation 
13  indicates  that,  for  a  reference  frame  moving  with  the  wind,  the  Bmnt-Vaisala  frequency  is  unaltered.  Flowever, 
measurements  are  made  in  a  stationary  reference  frame,  indicating  that  the  winds  control  the  measured  velocity  of 
the  gravity  waves.  The  buoyancy  frequency  depends  on  the  vertical  gradient  of  the  ambient  density.  From  the 
hydrostatic  equation  Po(z)=  -  p  og  and  the  ideal  gas  law  5po/5z  =  p  qR  T,  where  R  is  a  gas  constant  and  T  is  the 
temperature  in  degrees  Kelvin,  the  ambient  density  is  derived  as  p  o(z)  =  p  o(0)  exp(-g  z/RT)/RT.  The  vertical 
derivative  of  the  density  may  be  derived  from  this  relation. 

Equations  8,  10,  and  1 1  form  a  complete  set  of  first  order  linear  equations  for  the  pressure,  acoustic  particle  velocity 
and  density  fluctuations  associated  with  a  propagating  infrasound  wave  in  a  windy  medium,  valid  under  the 
assumptions  stated  previously.  To  summarize,  these  equations  are  valid  under  the  assumptions  that  the  ambient 
pressure,  density,  and  wind  fields  do  not  vary  with  time,  the  winds  are  horizontal,  and  vary  only  with  altitude.  That 
is,  the  wind  must  be  range-invariant. 

Finite  difference  modeling  of  infrasound  in  a  windy  environment 

Flere,  an  FDTD  method  is  outlined  for  infrasound  propagation  in  a  windy,  stratified  medium.  That  is,  it  is  assumed 
that  the  wind  blows  horizontally  and  is  invariant  along  range.  The  densities,  and  sound  velocities  vary  much  more 
gradually  in  the  horizontal  direction  than  vertically.  The  method  is  applied  to  a  2-D  model,  that  is,  the  source  is 
assumed  to  be  an  infinite  line  source.  The  method  resembles  an  FDTD  method  for  sound  waves  in  a  windy 
atmosphere  described  in  detail  by  Ostashev  et.al.  (2005),  but  with  the  addition  of  gravity. 


For  a  2-D  model  in  the  x-z  plane  with  horizontal  wind  speed  Wx,  Eqs.  (8)  becomes 

fhv  (h'x  ()wx 


Eq.  (10)  becomes 

and  Eq.  (11)  becomes 


f/l'i  OCx  OWx 

dt  “  dz 

dv~  dv^  -ifdp.^  \ 


dp„  dp„ 


dps 

dt 


=  ...  „ ,.  , ... 


dz 


dx  c'^  \  dt 


1  /op,,  (Jp.,\ 

^Kdt  Ox)' 


(14«) 

(Ub) 

(ir.) 

(10) 


The  static  sound  speed  c  and  ambient  density  p  o  may  vary  with  both  altitude  and  range;  the  wind  speed  Wx  varies 
only  with  altitude.  Equations.  14-16  are  in  a  form  suitable  to  computation  by  FDTD  techniques.  Note  that  for 
Wx  =  0,  that  is,  zero  wind  velocity,  and  g=0,  these  equations  reduce  to  the  usual  equations  for  acoustic  propagation  in 
a  static  medium  (e.g.,  Botteldoren,  1994). 

The  finite  difference  (FD)  method  relies  on  replacing  linear  partial  differential  equations  by  a  set  of  discrete 
equivalents.  Field  solutions  are  then  computed  over  a  discrete  set  of  nodes  that  comprise  the  spatial  grid. 

Figure  1  indicates  how  field  variables  are  defined  and  how  the  medium  is  discretized  for  the  staggered  grid  method, 
initially  developed  by  Yee  (1966).  The  model  is  decomposed  into  a  set  of  discrete  cells  of  dimension  Axx  Az.  The 
sound  speed  c  and  ambient  density  p  o  are  uniform  over  a  given  cell,  but  may  vary  from  cell  to  cell.  Pressure  and 
gravity  nodes  are  defined  at  the  center  of  each  cell  and  the  velocity  variables  are  located  midway  between  the 
pressure  nodes.  The  staggered  grid  formulation  increases  the  accuracy  of  the  FD  solution,  since  central  differences 
are  used  to  compute  the  discrete  derivatives  (Taflove  and  Flagness,  2000).  The  locations  of  the  vertical  velocity 
nodes  are  defined  in  such  a  way  as  to  allow  a  rigid  surface  (Vz=0)  to  be  defined  at  the  bottom  of  the  model. 
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Fignre  1.  The  finite  difference  model  is  decomposed  into  a  set  of  discrete  cells,  indicated  by  the  solid  lines, 
each  with  nniform  velocity  c  and  ambient  density  p  The  aconstic  pressnre  and  density 
pertnrbations  are  defined  by  the  nodes  at  the  center  of  each  cell,  indicated  by  the  filled  circles.  The 
locations  of  the  horizontal  velocity  Vr  (triangles  pointing  right)  and  vertical  velocity  variables  v^ 
(triangles  pointing  np)  are  defined  on  a  spatially  staggered  grid,  as  shown. 

Typically,  in  a  nonmoving  medium,  the  acoustic  velocities  and  pressures  are  computed  in  a  leap-frog  manner 
(Yee,  1966),  thus  the  velocities  and  pressures  are  computed  at  alternating  time-steps,  and  the  fields  from  the 
previous  time  step  are  overwritten.  Without  the  wind  terms,  the  governing  equations  (Equations  14-16)  indicate  that 
time  derivatives  in  the  acoustic  particle  velocities  depend  on  the  spatial  derivatives  of  the  pressure  variables,  and 
vice  versa.  However,  with  the  inclusion  of  advection  terms,  first  order  derivatives  in  space  and  time  must  be 
computed  simultaneously.  Various  numerical  implementations  have  been  suggested  for  the  computation  of  the 
FDTD  equations  for  sound  propagation  in  windy  environments  (Blumrich  and  Heinmann,  2002;  Van  Renterghem 
and  Botteldoren,  2003;  Ostashev  et.al.,  2005).  Here,  the  method  of  Ostashev  et.al.  (2005)  is  followed.  That  is, 
pressure  and  velocity  fields  are  over  saved  over  two  time  steps  so  that  time  derivatives  may  be  computed  using 
central  differences.  Refer  to  Ostashev  et.al.,  (2005)  for  further  detail. 

A  method  of  applying  FD  methods  in  the  frequency  domain,  in  cylindrical  coordinates,  for  a  point  source  in  an 
azimuthally  symmetric  model  was  developed  by  de  Groot-Hedlin  (2006). 

Comparison  ofFDFD  and  FDTD  solutions  in  a  model  with  topography 

In  this  subsection,  comparisons  are  made  between  FDFD  and  FDTD  results  for  a  shallow  underground  source  in  a 
medium  with  topography.  The  source  is  embedded  within  a  medium  with  a  sound  speed  2  km/s  and  a  density  of 
2000  k/m3.  The  sound  speed  within  the  air  is  shown  in  Figure  2. 
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Fignre  2.  Sonnd  velocity  profile  within  the  air.  This  sonnd  velocity  profile  is  taken  from  a  G2S  model  for  the 
Mt.  St.  Helens  region  for  March  9,  2005.  The  density  decreases  with  altitnde  hy  slightly  over  an 
order  of  magnitnde  over  this  altitnde  scale. 

The  sound  velocity  profile  shows  a  decrease  in  sound  speed  with  altitude  for  the  first  10  km;  this  corresponds  to  a 
typical  sound  velocity  profile  in  which  the  temperature  decreases  with  altitude.  The  effect  is  to  cause  some 
deflection  of  sound  upwards,  away  from  the  surface. 

The  finite  difference  time  domain  method  is  applied  in  Cartesian  coordinates,  thus  the  model  is  invariant  along  a 
direction  perpendicular  to  the  x-z  plane.  The  source  is  thus  a  line  source. 

This  model  features  a  broad,  symmetric  peak,  with  a  highest  altitude  of  2  km.  The  source  was  located  directly 
beneath  the  peak  at  an  altitude  of  1.2  km,  i.e.,  a  distance  of  0.9  km  below  the  ground  surface.  The  center  frequency 
of  the  source  was  0.5  Hz.  “Receivers”  were  located  at  intervals  of  5  km  from  each  side  of  the  peak  at  distances  from 
5  km  to  30  km  from  the  center  of  the  model. 

Several  “snapshots”  of  the  acoustic  pressure  are  shown  in  Figure  3.  As  indicated,  the  pressure  propagates  quickly 
through  the  ground  then  couples  to  the  air.  A  later  arrival  corresponds  to  acoustic  energy  that  couples  to  the  air  near 
the  source  and  propagates  outward  from  there.  Thus  there  should  be  two  main  arrivals  at  each  receiver,  a  first  one 
corresponding  to  acoustic  energy  propagating  through  the  earth,  coupling  to  the  air  near  the  receiver,  and  a  second 
corresponding  to  coupling  to  the  air  near  the  source.  The  traces  corresponding  to  this  model  are  shown  in  Figure  4. 
For  this  model,  ray -based  propagation  methods  which  rely  on  a  high  frequency  approximation  indicate  that  acoustic 
energy  is  refracted  upward,  away  from  the  ground  so  that  the  more  distant  receivers  lie  in  a  shadow  zone.  These 
results  thus  indicate  the  extent  to  which  synthesis  of  the  whole  waveform  is  required. 
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Fignre  3.  “Snapshots”  of  the  aconstic  pressnre  emanating  from  a  0.5  Hz  sonrce  at  0.9km  helow  the  central 
peak  at  the  time  points  5  s,  10  s,  and  15  s.  An  identical  color  scale  is  shown  for  each  plot. 
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diTiefs] 

Fignre  4.  Waveforms  at  distances  of  5  km  to  30  km  from  the  peak. 


range  [km] 

Fignre  5.  Transmission  loss  resnlts  derived  from  finite  difference  freqnency  domain  modeling  of  (top)  a 

0.5Hz  sonrce  at  0.6km,  i.e.,  1.5km  helow  the  snrface  and  (bottom)  for  the  same  sonrce  at  1.2km.  The 
second  sonrce  is  eqnivalent  to  the  one  shown  for  the  FDTD  simnlation.  Note  the  difference  in 
colorscales  for  the  two  sonrce  depths. 
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For  the  finite  difference  time  domain  results,  there  was  no  discemable  difference  in  the  results  with  or  without 
gravity,  as  gravitational  results  become  significant  at  much  lower  frequencies.  Incorporation  of  a  realistic  wind 
profile  showed  only  minor  asymmetry  in  the  waveform  results  on  either  side  of  the  peak.  Wind  and  gravity  have  not 
yet  been  incorporated  into  finite  difference  frequency  domain  modeling. 


CONCLUSIONS  AND  RECOMMENDATIONS 


Equations  have  been  developed  to  incorporate  the  effects  of  both  wind  and  gravity  into  an  FDTD  modeling  method. 
FDFD  methods  of  incorporating  these  effects  are  less  robust  as  they  require  the  inversion  of  large  matrices. 

It  has  been  demonstrated  that  a  full  waveform  modeling  technique  is  required  to  model  inffasound  generated  by  an 
underground  source.  Further  verification  of  the  FD  method  is  needed  for  sources  within  a  fraction  of  a  wavelength 
from  a  boundary  with  a  strong  impedance  contrast. 
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